NASA/TP-2010-216199 



An Improved Neutron Transport Algorithm for 
HZETRN 


Tony C. Slaba 

Old Dominion University, Norfolk, Virginia 

Steve R. Blattnig and Martha S. Clowdsley 
Langley Research Center, Hampton, Virginia 

Steven A. Walker 

Old Dominion University, Norfolk, Virginia 
Francis F. Badavi 

Christopher Newport University, Newport News, Virginia 


February 2010 




NASA STI Program ... in Profile 


Since its founding, NASA has been dedicated to 
the advancement of aeronautics and space science. 

The NASA scientific and technical information (STI) 
program plays a key part in helping NASA maintain 
this important role. 

The NASA STI program operates under the 
auspices of the Agency Chief Information Officer. It 
collects, organizes, provides for archiving, and 
disseminates NASA’s STI. The NASA STI program 
provides access to the NASA Aeronautics and Space 
Database and its public interface, the NASA Technical 
Report Server, thus providing one of the largest 
collections of aeronautical and space science STI in 
the world. Results are published in both non-NASA 
channels and by NASA in the NASA STI Report 
Series, which includes the following report types: 

• TECHNICAL PUBLICATION. Reports of 
completed research or a major significant phase 
of research that present the results of NASA 
programs and include extensive data or 
theoretical analysis. Includes compilations of 
significant scientific and technical data and 
infonnation deemed to be of continuing 
reference value. NASA counterpart of peer- 
reviewed formal professional papers, but having 
less stringent limitations on manuscript length 
and extent of graphic presentations. 

• TECHNICAL MEMORANDUM. Scientific 
and technical findings that are preliminary or of 
specialized interest, e.g., quick release reports, 
working papers, and bibliographies that contain 
minimal annotation. Does not contain extensive 
analysis. 

• CONTRACTOR REPORT. Scientific and 
technical findings by NASA-sponsored 
contractors and grantees. 


• CONFERENCE PUBLICATION. Collected 
papers from scientific and technical 
conferences, symposia, seminars, or other 
meetings sponsored or co-sponsored by NASA. 

• SPECIAL PUBLICATION. Scientific, 
technical, or historical information from NASA 
programs, projects, and missions, often 
concerned with subjects having substantial 
public interest. 

• TECHNICAL TRANSLATION. English- 
language translations of foreign scientific and 
technical material pertinent to NASA’s mission. 

Specialized services also include creating custom 

thesauri, building customized databases, and 

organizing and publishing research results. 

For more information about the NASA STI 

program, see the following: 

• Access the NASA STI program home page at 
http://www.sti. nasa.gov 

• E-mail your question via the Internet to 
help@sti.nasa.gov 

• Fax your question to the NASA STI Help Desk 
at 443-757-5803 

• Phone the NASA STI Help Desk at 
443-757-5802 

• Write to: 

NASA STI Help Desk 

NASA Center for AeroSpace Information 

7115 Standard Drive 

Hanover, MD 21076-1320 


NASA/TP-2010-216199 



An Improved Neutron Transport Algorithm for 
HZETRN 


Tony C. Slaba 

Old Dominion University, Norfolk, Virginia 

Steve R. Biattnig and Martha S. Clowdsley 
Langley Research Center, Hampton, Virginia 

Steven A. Walker 

Old Dominion University, Norfolk, Virginia 
Francis F. Badavi 

Christopher Newport University, Newport News, Virginia 


National Aeronautics and 
Space Administration 

Langley Research Center 
Hampton, Virginia 23681-2199 


February 2010 




Available from: 

NASA Center for AeroSpace Information 
7115 Standard Drive 
Hanover, MD 21076-1320 
443-757-5802 



Abbreviations 


AV 

Col 

GCR 

HZETRN 

SPE 


Average Value method 

Collocation Method 

Galactic Cosmic Rays 

High charge (Z) and Energy TRaN sport 

Solar Particle Events 



Contents 


Abstract 1 

1. Introduction 1 

2. Preliminary Energy Grid Convergence Testing 2 

3. Description of HZETRN Formalism 5 

4. Neutron Transport 8 

5. Results and Discussion 13 

6. Conclusions 22 

7. References 22 



Figures 


1 . 

2 . 

3. 

4. 

5. 

6 . 
7. 


9. 

10 . 
11 . 
12 . 

13. 

14. 

15. 

16. 

17. 

18. 

19. 

20 . 
21 . 


Dose versus depth in a 200 g/cm 2 Aluminum slab exposed to the August, 1972 King SPE spectrum 

using various energy grids 2 

Percent difference from average dose in a 50 g/cm 2 Aluminum slab exposed to the August, 1972 

King SPE spectrum using various energy grids 3 

Neutron fluence at 200 g/cnr in an Aluminum slab exposed to the August, 1972 King SPE spectrum 

using various energy grids 3 

Dose versus depth in a 200 g/cm 2 Aluminum slab exposed to the August, 1972 King SPE spectrum 

with elastic neutron production neglected in the transport procedure (All curves overlapping) 4 

Neutron fluence at 200 g/cnr in an Aluminum slab exposed to the August, 1972 King SPE spectrum 
with nuclear elastic neutron production neglected in the transport procedure (All curves 

overlapping) 4 

Comparison of scaled proton range ( R p (E)/vj ) and 4 He range ( Rj(E )) in Aluminum 6 

Neutron elastic production cross section for various projectile neutron energies and an Aluminum 

target 11 

Neutron fluence multiplied by proton stopping power and exponential attenuation term at 30 g/cm 2 

in an Aluminum target exposed to the August, 1972 King SPE spectrum 12 

Neutron fluence at 200 g/cm 2 in Aluminum exposed to the August, 1972 King SPE spectrum using 
various energy grids and two numerical methods for handling neutron elastic interactions (All curves 

overlapping) 13 

Neutron fluence at 200 g/cm 2 in water exposed to the August, 1972 King SPE spectrum using 
various energy grids and two numerical methods for handling neutron elastic interactions (All curves 

overlapping) 14 

Neutron flux at 200 g/cm 2 in Aluminum exposed to the 1977 solar minimum GCR spectrum using 
various energy grids and two numerical methods for handling neutron elastic interactions (All curves 

overlapping) 14 

Neutron flux at 200 g/cm 2 in water exposed to the 1977 solar minimum GCR spectrum using various 
energy grids and two numerical methods for handling neutron elastic interactions (All curves 

overlapping) 15 

Percent difference from converged neutron fluence at 200 g/cm 2 in Aluminum exposed to the 
August, 1972 King SPE spectrum using various energy grids and two numerical methods for 

handling neutron elastic interactions 16 

Percent difference from converged neutron fluence at 200 g/cm 2 in water exposed to the August, 

1972 King SPE spectrum using various energy grids and two numerical methods for handling 

neutron elastic interactions 16 

Percent difference from converged neutron flux at 200 g/cm 2 in Aluminum exposed to the 1977 solar 
minimum GCR spectrum using various energy grids and two numerical methods for handling 

neutron elastic interactions 17 

Percent difference from converged neutron flux at 200 g/cm 2 in water exposed to the 1977 solar 
minimum GCR spectrum using various energy grids and two numerical methods for handling 

neutron elastic interactions 17 

Dose versus depth in a 200 g/cm 2 Aluminum slab exposed to the August, 1972 King SPE spectrum 
using various energy grids and two numerical methods for handling neutron elastic interactions. ... 18 
Dose versus depth in a 200 g/cm 2 water slab exposed to the August, 1972 King SPE spectrum using 

various energy grids and two numerical methods for handling neutron elastic interactions 18 

Dose versus depth in a 200 g/cm 2 Aluminum slab exposed to the 1977 solar minimum GCR 
spectrum using various energy grids and two numerical methods for handling neutron elastic 

interactions 19 

Dose versus depth in a 200 g/cm 2 water slab exposed to the 1977 solar minimum GCR spectrum 
using various energy grids and two numerical methods for handling neutron elastic interactions. ... 19 
Percent difference from converged dose spectrum in Aluminum exposed to the August, 1972 King 
SPE spectrum using various energy grids and two numerical methods for handling neutron elastic 
interactions 20 


v 



22. Percent difference from converged dose spectrum in water exposed to the August, 1972 King SPE 

spectrum using various energy grids and two numerical methods for handling neutron elastic 
interactions 20 

23. Percent difference from converged dose spectrum in Aluminum exposed to the 1977 solar minimum 

GCR spectrum using various energy grids and two numerical methods for handling neutron elastic 
interactions 21 

24. Percent difference from converged dose spectrum in water exposed to the 1977 solar minimum GCR 

spectrum using various energy grids and two numerical methods for handling neutron elastic 
interactions 21 


vi 



Tables 


1 . 


Energy domains of the neutron elastic production cross section for Hydrogen, Oxygen, and 
Aluminum targets for projectile kinetic energies of 0.01, 0.1, 10, and 100 MeV 10 



viii 



Abstract 


Long term human presence in space requires the inclusion of radiation constraints in 
mission planning and the design of shielding materials, structures, and vehicles. It is 
therefore necessary to expose the numerical tools commonly used in radiation analyses to 
extensive verification, validation, and uncertainty quantification. In this paper, the 
numerical error associated with energy discretization in HZETRN is addressed. An 
inadequate numerical integration scheme in the transport algorithm is shown to produce 
large errors in the low energy’ portion of the neutron and light ion fluence spectra. It is 
further shown that the errors result from the narrow energy domain of the neutron elastic 
cross section spectral distributions, and that an extremely fine energy grid is required to 
resolve the problem under the current formulation. Since adding a sufficient number of 
energy’ points will render the code computationally inefficient, we revisit the light ion and 
neutron transport theory developed for HZETRN and focus on neutron elastic 
interactions. Two numerical methods are developed to provide adequate resolution in the 
energy domain and more accurately resolve the neutron elastic interactions. 
Convergence testing is completed by running the code for various environments and 
shielding materials with various energy’ grids to ensure stability of the newly 
implemented method. 

1. Introduction 

In free space and low Earth orbit, it is well known that the dose within the first few centimeters of 
shielding in a space vehicle, space station, or habitat is dominated by energy deposition from high energy 
charged particles, while at larger depths the dose is increasingly dependent on secondary neutron 
production [Getselev et al. 2004]. This implies that for space structures, where the cumulative shielding 
thicknesses can be substantial, accurate dose estimates depend heavily on accurate descriptions of the 
neutron fluence spectrum. Since little neutron data exists for these environments and shielding 
configurations, code validation is particularly difficult. However, verification and uncertainty quantification 
are still necessary, and convergence testing is a necessary step in both efforts. 

The deterministic transport code HZETRN (High charge (Z) and Energy TRaNsport) [Wilson et 
al. 1986, 1989, 1991, 2006; Cucinotta et al. 1993; Shinn et al. 1991], developed at NASA Langley 
Research Center, obtains an approximate solution to the Boltzmann transport equation within the straight 
ahead and continuous slowing down approximations. While the code has been documented in its accuracy 
and applicability as a design tool in assessing space radiation exposure [Wilson et al. 2005], convergence 
testing of the light ion (A < 4, where A is the atomic mass of the particle) and neutron transport routines has 
been limited to evaluation of dosimetric quantities using three energy grids in various materials exposed to 
solar particle event (SPE) environments [Shinn et al. 1991]. Such testing is certainly a necessary step in 
evaluating code accuracy but does not guarantee convergence of individual ion fluence. Further, interest in 
a fluence based approach to risk assessment suggests that dose and dose equivalent may not be sufficient in 
evaluating code accuracy and that individual ion fluence should be used instead [Cucinotta et al. 2006]. 

In this work, we first show that large errors are introduced in dose estimates past 50 g/cm 2 of 
shielding by altering the character and fidelity of the energy grid used. The errors are determined to be 
related to the neutron elastic production cross section and are highly visible in the low energy neutron 
fluence spectrum. These errors had not been seen in previous studies [Shinn et al. 1991] due to the low 
fidelity of the energy grids studied. One could resolve the problem by sufficiently increasing the fidelity of 
the energy grid, but such a refinement would render the code computationally inefficient. Therefore, 
alternative methods need to be considered. In order to develop the most appropriate solution to this 
problem, we review the light ion and neutron transport procedure developed for HZETRN [Wilson et al. 
1989, 1991, 2006; Cucinotta et al. 1993; Shinn et al. 1991] and give special attention to low energy 
neutrons. Two numerical procedures are presented, and they are shown to produce stable low energy 
neutron spectra with fewer energy grid points than what was used previously. 
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2. Preliminary Energy Grid Convergence Testing 


In this section, we examine energy discretization error in HZETRN by evaluating the code with 
several different energy grids. A 200 g/cm 2 aluminum slab was exposed to the August, 1972 King [King et 
al. 1974] SPE spectrum. Twelve different energy grids were used for the transport procedure. The first six 
were equally log spaced in energy and contained 50, 100, 150, 200, 300, and 500 points (E-50, E-100, E- 
150, E-200, E-300, E-500 in the plots below). The remaining six were equally log spaced in proton range 
and also contained 50, 100, 150, 200, 300, and 500 points (R-50, R-100, R-150, R-200, R-300, R-500 in the 
plots below). It will be shown later that the transport procedures in HZETRN are defined in terms of proton 
range, and so testing both of these classes allowed for a more detailed comparison. Note also that the 
conversion from proton range to energy is accomplished using the range-energy relation 


r 


r E dE 

Jo S(E')’ 


( 1 ) 


where r(E) and S(E) are the proton range and stopping power, respectively, and E is the proton kinetic 
energy. 

In Figure 1, it is clear that large differences are introduced into dose values for depths larger than 
50 g/cm 2 by using different energy grids in the transport procedure. Since the vertical axis of Figure 1 
covers six orders of magnitude in log scale, Figure 2 is given to more clearly quantify the differences at 
depths less than 50 g/cm 2 . The average of all the dose estimates generated by each energy grid was 
calculated at discrete depths, and Figure 2 gives the percent difference of each dose estimate from this 
average. 



Figure 1. Dose versus depth in a 200 g/cm 2 Aluminum slab exposed to the August, 1972 King SPE 

spectrum using various energy grids. 
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2 

Depth in Aluminum (g/cnT) 


Figure 2. Percent difference from average dose in a 50 g/cm 2 Aluminum slab exposed to the August, 1972 

King SPE spectrum using various energy grids. 



Energy (MeV) 


Figure 3. Neutron fluence at 200 g/cm' in an Aluminum slab exposed to the August, 1972 King SPE 

spectrum using various energy grids. 
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Figure 4. Dose versus depth in a 200 g/cm 2 Aluminum slab exposed to the August, 1972 King SPE 
spectrum with elastic neutron production neglected in the transport procedure (All curves overlapping). 



Energy (MeV) 


Figure 5. Neutron fluence at 200 g/cm" in an Aluminum slab exposed to the August, 1972 King SPE 
spectrum with nuclear elastic neutron production neglected in the transport procedure (All curves 

overlapping). 
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In order to find the source of these errors, it was first noted that for an SPE environment, 
HZETRN transports neutrons, protons, 2 H, 3 H, 'He, and 4 He. These fluence spectra were examined in 
detail, and it was found that the largest errors appeared in the low energy ( < 30 MeV) neutron spectrum, as 
seen in Figure 3. The remaining light ions also displayed some sensitivity to the choice of energy grid, but 
this was hypothesized to be a result of the large fluctuations seen in the neutrons. 

There are several neutron production channels included in HZETRN, and elastic interactions are 
the dominant ones for neutron energies less than 10 MeV [Wilson et al. 1991]. This can be verified 
numerically by neglecting such processes in the transport procedure. Figures 4 and 5 give the dose as a 
function of depth in the Aluminum slab and the neutron spectrum at 200 g/cnr in an Aluminum slab 
exposed to the August, 1972 King SPE spectrum with elastic neutron production terms being neglected in 
the transport procedure. Figure 4 and Figure 5 illustrate that neglecting the elastic neutron production terms 
in the transport procedure removed most of the instabilities seen in Figure 1 and Figure 2. In Figure 5, the 
neutron spectra are so close that one cannot discern among them. It is concluded that either the elastic 
neutron production terms are highly sensitive to the energy grid used, or the numerical transport procedure 
in HZETRN is not adequately handling such interactions. We now proceed to revisit the neutron and light 
ion transport procedure used in HZETRN. 

3. Description of H ZE TRN Formalism 

In this section, we review the formulation of the numerical marching procedures used in 
HZETRN. The Boltzmann transport equation with the continuous slowing down and straight ahead 
approximations is given as [Wilson et al. 1991] 



with the linear differential operator 


m = 


d_ 

dx 


J__d_ 
A. dE 

j 


Sj(E) + a ■ (E) 


(x,E), 


(2) 


( 3 ) 


where <Pj(x,E) is the fluence of type j particles at depth x with kinetic energy E, Aj is the atomic mass 
number of a type j particle, Sj(E) is the stopping power of a type j ion with kinetic energy E, g{E) is the 
total macroscopic cross section for a type j particle with kinetic energy E, and ajAE,E') is the macroscopic 
production cross section for interactions in which a type k particle with kinetic energy E' produce a type j 
particle with kinetic energy E. 

The continuous slowing down operator in equation (3) 


J_ _d_ 
A] dE 



( 4 ) 


is an approximate representation of the rate at which charged particles lose energy due to atomic 
interactions with target material. Though atomic interactions cause charged particles to lose energy in 
discrete increments as they pass through a material, there are sufficiently many of these interactions in a 
unit path length to justify a continuous approximation [Wilson et al. 1991]. It is advantageous to further 
approximate this term by considering the relation from Bethe stopping power theory [Bethe et al. 1930] 


v r 
j j 


( 5 ) 
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Energy (MeV/amu) 


Figure 6. Comparison of scaled proton range ( R p (E)/vj ) and 4 He range (Rj(E)) in Aluminum. 


where u- = Z 1 - / A • , Z, is the atomic charge of a type j particle, and is defined by the range-energy 
relation 


r.(E) = Aj‘ 


E dE' 


S 3 {E' 


(6) 


Equations (5) and (6) suggest that 


u^E) 


A j 


( 7 ) 


where S(E) is the proton stopping power. The approximation in equation (5) is less accurate at low 
energies, or low residual range, as seen in Figure 6, but it is expected to have a small effect on final dose 
estimates [Wilson et al. 2006]. 

Substitution of equation (7) into the continuous slowing down operator implies that the linear 
differential operator given by equation (3) is now written 


m 


d__ d_ 
dx V ' J dE 


S(E ) + aAE) 


so that upon introduction of the scaled quantities 


i’ji.x, r) = v -S(E)<f)-(x, E ) , 


( 8 ) 


( 9 ) 


s jk ( r > r ') = S{E)a - k (E, E ') , 


( 10 ) 
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the transport equation can be represented in terms of the single range variable 


d d 

V b <7 

dx 3 dr 3 


(r) 


IS ■ n oo 

i’j (x, r) = 2_^ ~ J r s jk ( r , r ') V’fc {x, r ') dr 1 , 


( 11 ) 


where r is the residual proton range given in equation (1). It is noted for neutrons, that v- is taken to be 

unity in fluence scaling relations and formally taken as zero in range scaling relations; more will be said 
about this shortly. 

Equation (11) can be inverted using the method of characteristics [Haberman 1998] to obtain the 
Volterra type integral equation [Wilson et al. 2006] 

ipj (x, r) = ,r + t'jx) 

_ IS ■ r*x r*oo _ c ( r „t\ (12) 

+Y^ — J 0 J , e 3 ’ s jk( r + u j x \r')i> k {x - x\r') dr'dx' , 
k ^ k r V i X 

with the exponential term 

Sj(r,x) = f o a^r + vfidt , (13) 

Note that Vj in equation (9) and Vj/ v k in equation (11) are both fluence scaling relations, and v„ = 1 in both 
cases to provide a nontrivial scaling. Conversely, wherever v appears as the argument of a fluence or cross 
section in equations (12) and (13), it is a range scaling relation and v„ = 0. This convention is taken to 
reflect the absence of atomic interactions in neutron transport ( S„ (E) = 0 ). 

Marching procedures are developed for equation (12) by considering light ions (A < 4) and heavy 
ions (A > 4) separately as demanded by the nature of the nuclear cross sections and boundary conditions. 
For heavy ions, it is noted that projectile fragments have energy and direction very near that of the 
projectile, while target fragments are produced nearly isotropically with very low energy [Wilson et al. 
2005]. Though projectile and target fragments cannot be explicitly decomposed in this manner, it has been 
shown that an approximate decoupling is suitable for space applications [Wilson et al. 1991]. These 
observations allow heavy target fragments to be neglected in the transport procedure [Wilson et al. 1989, 
1991] (their contribution to dose is approximately accounted for after the transport procedure) and the 
scaled production cross section appearing in equation (12) to be recast as 

s jk (r,r') = a jk (r)6(r - r') , (14) 

to reflect the nearly equal energies (range) of the projectile and fragment. The final marching procedure for 
heavy ions was developed by Wilson et al. [1986] 


ipjix + h,r) = e ^ r ’ h ^'ip j (x,r + v,.h) + ^2 — cj 


k>j ^ k 


jk 


r + is 


~ a > r+I/ 4r ” CT t( r +(^+ 1 't)y ft 
e K — e K ZJ 


r + i v j + u k) g 


!N 

0{h 2 ), 


x,r + {v , +f k )- 


(15) 


r + v- 


where h is the marching step-size in g/cm 2 . The summation over k > j reflects the fact that projectile 
fragments will always have mass less than the original projectile. 
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Alternatively, for light ions, both projectile and target fragments are included in the transport 
procedure, and the broad energy distribution in collision events indicates that equation (14) cannot be used. 
The final light ion marching procedure was also developed by Wilson et al. [2006] 


'ipjix + h,r) = e ^ r,h) 'ip j (x,r + z^/i) + J2~f l 




k v h J r+u^i% 


xF£(r,r';h)ip k 


x, r '+ v, - 

k 2 


dr'+ 0(h 2 ), 


(16) 


where h is the marching step-size in g/cm 2 , and the summation over k is taken over all light ions. Note that 
the integrand has been simplified using 


Ffc(r, r '-h) = f S j k (r + u jZ ,r')dz. 


(17) 


The nature of the boundary condition determines how equations (15) and (16) are evaluated and 
coupled. For most SPE and trapped proton environments, there is a negligible heavy ion component, and 
only the light ion marching procedure (equation (16)) is evaluated. For galactic cosmic ray (GCR) 
environments, the heavy ion component is transported using equation (15), and the light ion component is 
transported using a modified form of equation (16), namely 


tpAx + h,r) = e iil - r ’ h ' > ip j (x,r + v.h) 


-E- r 

k V k Jr +v,h/2 


e HSl C T ’ 2 J F£(r,r';h)ip k 


-EE f 


k>J v k 


a jk I r + V ; ~ 


A 


x,r + {v ■ + v k ) — 


x, r '+ u k - 


dr' 


~ a J r + ^y _ ~ a i.\r+^ j +i' t ) l ^ h 


r +Ai+ V k) 


h' 


— cr ■ 


r + v ■ 


0(h 2 ) 


(18) 


where the first summation over k is again taken over all light ions, and the second summation over k is 
taken over heavy ion projectiles where the lower limit of the summation, J, represents the heaviest light 
ion. The second summation physically accounts for light ion production from heavy ion projectiles. In both 
the light particle and heavy ion marching equations, the 0(h 2 ) terms are all small. Flereafter, they will not 
be written. 


4. Neutron Transport 

It was shown in the preliminary convergence testing that instabilities in the neutron transport 
procedure were caused by the elastic component of the neutron production cross section. This process is 

accounted for in the integrated light particle production cross sections F^(r,r';ti) (equation (17)), and so 

we must examine the light ion marching procedure to resolve this issue. 

With a simple shift of the integration variable r\ one can rewrite the light ion marching procedure, 
equation (16), as 



( 19 ) 


ipj(x + h,r) = e + vfr) 

+ -y *J_ r°° e -( J (r,h/2)-C k (r'+v J h/2,h/2) 
k "k Jr 

xF^ (r, r '+ Vjh / 2 ;/i) ip k (x,r'+ (y j +v k )h/2)dr' . 

Equation (19) is evaluated by first considering the following approximation from Wilson et al. [1991] for a 
fixed r m 

OO 

f r K ( r m y)^j( x ,r')dr' ttJ2 K (r m ,r l )f r ,+1 ip j ( x y) dr '’ (20) 

m l—m 1 


where r ; is the average value of r ; and r (+1 , and oo denotes some upper limit dependent on the energy 

spectrum of the boundary condition. The accuracy of this approximation depends both on the spacing of 
the r-grid as well as the behavior of the kernel over the region [r/,r /+1 ] for each / and r,„. Such an 
approximation requires that the kernel is nearly constant or slowly varying over the region [r/,r/+ j]. In order 
to apply equation (20), the attenuation terms are approximated, and for a fixed r m the light ion marching 
procedure is now evaluated as 


i'Ax + h, r ) = e a ’ lr " )h 'J’ i {x,r + I'M 


k l—m ^k 


"j V ’ m 
j ~ a j( r m '> h / 2 - a k( r l ) h /' 2 pA 


F £( r m ’ r l +V.h/2-h) 


x [$*(*,»} + Wj + V k )h/ 2) - $ k (x,r l+1 + (y 3 + v k )h/ 2)] , 

where we have used the following definition for integral fluence 

X oo noo 

4> j (x,r')dr' = u jJ E 4>j (x,E')dE' . 


( 21 ) 


( 22 ) 


For neutrons, v n = 0 in range scaling relations, and equation (21) takes the special form 

1> n ( x + ft ’ r J = e ~ an(rm)h ^M' r m) 

OO 

k^n l—m k 

x [$*(*>»/ + v k h / 2 ) - ®k( x ’ r i+i + u k h / 2 )] ( 23 ) 

OO 

l—m 

x [ $ n (*. r /)- $ n (*. r i+i)] • 


The elastic component of the neutron production cross section is contained in the term (r m , i ] ; h) , and 
so we explicitly decompose this term for neutron production from neutron-nucleus interactions 

F nn(Vvh) = S(e(rJ)[a™ ( e (rJ,e(ij)) + (e(rj,e(fj))] (24) 

where the superscript re denotes the nuclear reactive component, the superscript el denotes the nuclear 
elastic component, and e{r) is the energy associated with the proton range r. 
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Table 1. Energy domains of the neutron elastic production cross section for Hydrogen, Oxygen, and 
Aluminum targets for projectile kinetic energies of 0.01, 0.1, 10, and 100 MeV. 


a t 

a 

Scattered neutron energy ranges for fixed projectile energies (MeV) 

0.01 

0.1 

10 

100 

1 

0 

(0,0.01) 

(0,0.1) 

(0,10) 

(0,100) 

16 

0.7785 

(0.0079,0.01) 

(0.07785,0.1) 

(7.785,10) 

(77.85,100) 

27 

0.8622 

(0.00862,0.01) 

(0.0862,0.1) 

(8.62,10) 

(88.62,100) 


The neutron elastic production cross section is restricted by kinematics to be nonzero only for 
scattered neutron energies bounded by [Bell et al. 1970] 

E'a<E<E', (25) 


where the parameter a is [Bell et al. 1970] 


a = 


' A - if 

^ Aj, + 1 


(26) 


and A t is the atomic weight of the struck nucleus. Table 1 gives the kinematic regions for 0.01, 0.1, 10, and 
100 MeV neutron projectiles on Hydrogen, Oxygen, and Aluminum targets, respectively. While the regions 
are quite broad for Hydrogen targets, it is clear that heavier targets cause the regions to narrow 
considerably. 

From equations (20) and (23), it is clear that the following approximation is being made for 
neutron elastic interactions 


f e <T « (r ' )ft/2 CT ^ (£:(/■ ),e(r')) ip (x, r’)dr’ 

*J T m 

oo 


l—m 


(27) 


This equation can also be expressed completely in terms of energy as 



e -<rn( E ') h / 2 (j el 
u nn 


(E m ,E')S(E')<P n (x,E')dE 


I 


« Y. e ~ an{El)hl2 <n(KM*nMri)) - * n (x,e(r l+1 ))} 

l—m 


(28) 


where £) = (£j + E l+1 ) / 2 . We now examine the special case of an Aluminum target and E m = 0.01 

MeV. To use the approximation in equations (27) or (28), the energy grid must be fine enough so that there 
are at least two grid points in the region (0.01 - 0.116 MeV). For equal log spacing in energy, this requires 
at least 170 grid points - 70 points more than what has been previously used in HZETRN, and sufficiently 
many to render the code computationally inefficient. 
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Figure 7. Neutron elastic production cross section for various projectile neutron energies and an Aluminum 

target. 


Two alternative numerical approximations are now considered to resolve this problem. The 
average value approximation is very similar to the approximation in equation (28), and therefore, it fits 
easily into the existing transport procedure. It can be written as 


/” {E ' )h ' 2 °L (■ E m ,E')S(E') <f> n (x, E')dE' 

oo 

~ Z) e ~ an{El)h/2 - ® n M r i+ 1))] » 


(29) 


with the average value of the neutron elastic production cross section 


,(3 


m 5 E l 


E l + 1 E l 


f E (T el (E ,E')dE' 

J £ ran V m ’ / 


(30) 


evaluated using a 10 point Gauss-Legendre numerical integration scheme. The upper limit of integration 
E = min | E [+: . E m / a} is used to minimize the number of quadrature points. This approximation will 

be exact if the neutron production cross section behaves linearly over the region of integration. Figure 7 
shows the neutron elastic production cross section generated by F1ZETRN for various projectile neutron 
energies on an Aluminum target. It is clear that the cross section is nearly linear for projectile energies less 
than 0.1 MeV, but it is considerably more peaked for energies greater than 1 MeV. The average value 
approximation will not be as accurate in these regions, but it will be shown later that the errors are small in 
most cases. 

The collocation approximation is included here as a verification of the average value 
approximation. The integral in equation (28) is estimated according to 
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with the coefficients 


= j E E M BAE')cj* n {E m ,E')dE' , 


defined in terms of the linear basis splines 


B,(E) = 


(E ~ Ef-j) / (E t — £ G 

(-®t+i ~ E) / (E l+1 — E { ), E £ [£' ; ,£ l i+1 ] . 
0 , otherwise 


The lower and upper limits of integration 

E(-) = ma x{E l _ v E m }, 
E (+ ) = min {E l+1 ,E m / a), 


(32) 


(33) 


(34) 

(35) 


are chosen to ensure that the integration variable is simultaneously within both the kinematic region of the 
production cross section and the non-zero region of the basis spline. The coefficients a hn are also evaluated 
using a 10 point Gauss-Legendre numerical integration scheme. 



Energy (MeV) 


Figure 8. Neutron fluence multiplied by proton stopping power and exponential attenuation term at 30 
g/cm 2 in an Aluminum target exposed to the August, 1972 King SPE spectrum. 
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The collocation approximation becomes exact if the product of the exponential term, proton 
stopping power, and neutron fluence behave linearly over the regions of integration. Figure 8 shows this 
product at 30 g/cnr in an Aluminum target exposed to the August, 1972 King SPE spectrum. The scaled 
neutron fluence covers eight orders of magnitude over the specified energy domain, and so it is certainly 
not linear. However, this quantity will be integrated over very small energy regions over which the scaled 
neutron fluence varies slightly. We therefore expect the collocation approximation to be more accurate than 
the average value approximation. 

In the following section, the accuracy of the proposed numerical solutions will be tested by 
repeating the energy grid convergence study presented earlier for large Aluminum and water targets 
exposed to both an SPE and GCR environment. 

5. Results and Discussion 

In this section, we present the results of convergence studies comparing the two new numerical 
methods for neutron transport. As before, two classes of energy grids were used for the transport 
procedures (logarithmic spacing in energy and logarithmic spacing in proton range) in combination with 
three different sizes (100, 300, and 500 points) and the two different numerical procedures for handling 
neutron elastic interactions outlined in the previous section. In the plots below, the letters E (logarithmic 
spacing in energy) and R (logarithmic spacing in proton range) have been used to identify the classes of 
grids used, the numerical value (100, 300, or 500) identifies the number of grid points used, and (AV) or 
(Col) identifies the method used to handle neutron elastic interactions (average value or collocation, 
respectively). 

Figures 9-12 give the neutron fluence spectra at 200 g/cnr in Aluminum and water targets 
exposed to the August, 1972 King SPE and 1977 solar minimum [O’Neill et al. 2006] GCR spectra, 
respectively. 



Energy (MeV) 


Figure 9. Neutron fluence at 200 g/cm 2 in Aluminum exposed to the August, 1972 King SPE spectrum 
using various energy grids and two numerical methods for handling neutron elastic interactions (All curves 

overlapping). 
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Energy (MeV) 


Figure 10. Neutron fluence at 200 g/cm 2 in water exposed to the August, 1972 King SPE spectrum using 
various energy grids and two numerical methods for handling neutron elastic interactions (All curves 

overlapping). 



Energy (MeV) 


Figure 11. Neutron flux at 200 g/cm 2 in Aluminum exposed to the 1977 solar minimum GCR spectrum 
using various energy grids and two numerical methods for handling neutron elastic interactions (All curves 

overlapping). 
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Energy (MeV) 

Figure 12. Neutron flux at 200 g/cm 2 in water exposed to the 1977 solar minimum GCR spectrum using 
various energy grids and two numerical methods for handling neutron elastic interactions (All curves 

overlapping). 


Figures 9-12 show that the large errors seen previously in the low energy portion of the neutron 
fluence spectra have been resolved. The proposed numerical methods agree very well at the lowest energies 
and at large depths. The remaining error can be quantified by identifying a converged solution and 
comparing all other results to the converged result. We define the converged results as those obtained with 
the collocation method using the R-500 grid. This choice was motivated by studying the data represented in 
Figures 9-12 and recognizing that the collocation method converges more rapidly, in most cases, as a 
function of energy grid size than does the average value method. It is worth noting that the average value 
method with the R-500 grid agrees almost identically with the converged results. The percent difference of 
each neutron fluence (or flux) spectrum from the converged neutron fluence (or flux) spectrum was 
calculated, and the results are given in Figures 13-16. 

Note that in each of the Figures 13-16 there is significant error near the end of the energy 
spectrum. This error is caused by truncating the energy domain of the external space radiation enviromnent 
to a finite subset of its actual semi-infinite energy domain. The error is most evident in the last few grid 
points as their distance from the upper bound of the energy domain will change significantly if the fidelity 
of the energy grid is altered. It is important to note that neutron fluences at these energies are extremely 
small, and so these errors are unimportant for dose calculations. Looking at the remainder of the energy 
spectrum {E < 100 MeV), it appears that errors are bounded by 40% for the SPE environment and 25% for 
the GCR environment. In general, it appears that the collocation method is more accurate in Aluminum, 
while the average value method is slightly more accurate in water. 

Figures 17-20 give dose as a function of depth in Aluminum and water targets exposed to the 
August, 1972 King SPE spectrum and the 1977 solar minimum GCR spectrum, respectively. Figures 17-20 
show that the large errors introduced in dose estimates at large depths have also been greatly reduced, but 
non-negligible differences still exist at depths greater than 50 g/cm 2 . Note that the vertical axes of Figures 
17 and 18 are in log scale to more accurately depict the rapid dose attenuation common to SPE 
environments, while the vertical axes of Figures 19 and 20 are in linear scale to more accurately depict the 
moderate dose attenuation common to GCR environments. 
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Figure 13. Percent difference from converged neutron fluence at 200 g/cnr in Aluminum exposed to the 
August, 1972 King SPE spectrum using various energy grids and two numerical methods for handling 

neutron elastic interactions. 



Energy (MeV) 


Figure 14. Percent difference from converged neutron fluence at 200 g/cm 2 in water exposed to the August, 
1972 King SPE spectrum using various energy grids and two numerical methods for handling neutron 

elastic interactions. 
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Energy (MeV) 


Figure 15. Percent difference from converged neutron flux at 200 g/cm 2 in Aluminum exposed to the 1977 
solar minimum GCR spectrum using various energy grids and two numerical methods for handling neutron 

elastic interactions. 



Energy (MeV) 


Figure 16. Percent difference from converged neutron flux at 200 g/cm 2 in water exposed to the 1977 solar 
minimum GCR spectrum using various energy grids and two numerical methods for handling neutron 

elastic interactions. 
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Figure 17. Dose versus depth in a 200 g/cm Aluminum slab exposed to the August, 1972 King SPE 
spectrum using various energy grids and two numerical methods for handling neutron elastic interactions. 



Figure 18. Dose versus depth in a 200 g/cm 2 water slab exposed to the August, 1972 King SPE spectrum 
using various energy grids and two numerical methods for handling neutron elastic interactions. 
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Figure 19. Dose versus depth in a 200 g/cm 2 Aluminum slab exposed to the 1977 solar minimum GCR 
spectrum using various energy grids and two numerical methods for handling neutron elastic interactions. 



Figure 20. Dose versus depth in a 200 g/cm 2 water slab exposed to the 1977 solar minimum GCR spectrum 
using various energy grids and two numerical methods for handling neutron elastic interactions. 
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To quantify the remaining error, we again compare all results to the converged result. The percent 
difference of each dose spectrum from the converged dose spectrum was calculated, and the results are 
given in Figures 21-24. Note that in each of the Figures 21-24, the increasing error as a function of depth 
reflects the accumulated error due to energy discretization and is expected for the iterative transport 
algorithm used in HZETRN. The errors are bounded by 25% for the SPE environment and 15% for the 
GCR environment. 



Figure 21. Percent difference from converged dose spectrum in Aluminum exposed to the August, 1972 
King SPE spectrum using various energy grids and two numerical methods for handling neutron elastic 

interactions. 
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Figure 22. Percent difference from converged dose spectrum in water exposed to the August, 1972 King 
SPE spectrum using various energy grids and two numerical methods for handling neutron elastic 

interactions. 
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Figure 23. Percent difference from converged dose spectrum in Aluminum exposed to the 1977 solar 
minimum GCR spectrum using various energy grids and two numerical methods for handling neutron 

elastic interactions. 



Figure 24. Percent difference from converged dose spectrum in water exposed to the 1977 solar minimum 
GCR spectrum using various energy grids and two numerical methods for handling neutron elastic 

interactions. 
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6. Conclusions 


Energy grid convergence testing with HZETRN showed that signficant errors are introduced into 
low energy neutron fluence spectra and dose estimates for an Aluminum target exposed to the August, 1972 
King SPE spectrum. The errors were shown to be caused by a numerical procedure that was incapable of 
resolving the fine energy domain of the neutron elastic production cross sections. Two numerical methods 
were developed specifically to handle this problem and were shown to resolve the errors seen in the 
preliminary convergence testing. The methods were found to be comparable in most cases, and energy grid 
discretization error was shown to decrease properly as a function of increasing grid size. Though the largest 
errors were obtained using the 100 point energy grids, the errors on dose were bounded by 25% in SPE 
environments and 15% in GCR environments. Due to the increased computational efficiency achieved by 
the smaller grid structures, and the ease with which the average value method can be implemented, it is 
concluded that the average value method with 100 energy grid points is sufficient for future use in 
HZETRN. Future work should focus on reducing the remaining sources of energy grid discretization error, 
and consider coupled convergence studies in both step-size and energy grid size. 
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